TODO:

Formulas

Calculate the upwind area aka footprint (\(X_f\), m) representative of these measurements that accounts for a percentage of the flux. Where \(U\) is the assumed constant wind speed, \(u^{*}\) is the friction velocity, \(d\) is the displacement height due to vegetation (taken as 2/3 the vegetation height), \(z\) is the instrument height, \(k\) is the von karmans constant (usually taken as 0.4), \(f\) is the fraction of the area that explains the flux [@Schuepp1990]:

\[ X_f = - \frac{1}{ln(f)}\frac{U}{u^*}\frac{z-d}{k} \] Friction velocity (\(u_*\)) can be modelled as in wind_profile.Rmd or is also measured at the Eddy Covariance Sites.

find the parameters \(u^{∗}\), \(Z_{0m}\), and \(d_0\) of the logarithmic wind profile equation that best matches the measured data provided below:

Constants:

# need to measure this in the field still (2022-06-01)
temp_offset <- 125

k <- 0.4 # von karmans constant (unitless) 

d0 <- 1.965869 # m estimated from jan 16 2022 at 4:00 pm 
f <- 0.8 # 80%

Load Data:

low_wind <- readRDS('data/eddy_cov_15min_combined.rds') |> select(TIMESTAMP, wnd_dir_compass, u_star, low_ec_rslt_wnd_spd = result_wind_spd)# |> 
#   pivot_longer(u_star:low_ec_rslt_wnd_spd)
# 
# ggplot(low_wind, aes(TIMESTAMP, value, colour = name)) + geom_line()
# 
# plotly::ggplotly()

mid_hi_wind <- readRDS('../../analysis/treefort-main/data/met/met_main.rds') |> select(TIMESTAMP, USWindSpeed_S_WVT, top_ec_rlst_wnd_spd, top_ec_u_star)

wind_wide <- left_join(low_wind, mid_hi_wind, by = 'TIMESTAMP')

wind_dir <- wind_wide |> 
  mutate(wind_dir = (wnd_dir_compass + temp_offset) %% 365) |> 
  select(TIMESTAMP, wind_dir) 

low <- wind_wide |> 
  select(TIMESTAMP, u_star, wind_speed = low_ec_rslt_wnd_spd) |> 
  mutate(inst_height = 3) |> 
  mutate(
    med = zoo::rollapply(wind_speed, width = 1000, FUN = median, na.rm = T, align = 'center', partial = T),
    sd = zoo::rollapply(wind_speed, width = 1000, FUN = sd, na.rm = T, align = 'center', partial = T),
    stdep = (wind_speed - med) / sd,
    flag = ifelse(wind_speed > 4, T, F)
  ) |> 
    select(TIMESTAMP, u_star, wind_speed, inst_height, flag) 


ggplot(low, aes(TIMESTAMP, wind_speed, colour = flag)) + geom_point()
## Warning: Removed 1372 rows containing missing values (geom_point).

plotly::ggplotly()
hi <- wind_wide |> 
  select(TIMESTAMP, u_star = top_ec_u_star, wind_speed = top_ec_rlst_wnd_spd) |> 
  mutate(
    inst_height = 13.5,
    mean = zoo::rollapply(wind_speed, width = 48, FUN = mean, na.rm = T, align = 'center', partial = T),
    sd = zoo::rollapply(wind_speed, width = 48, FUN = sd, na.rm = T, align = 'center', partial = T),
    stdep = (wind_speed - mean) / sd,
    flag = ifelse(wind_speed > 20, T, F)
  ) |> 
  select(TIMESTAMP, u_star, wind_speed, inst_height, flag) 

ggplot(hi, aes(TIMESTAMP, wind_speed, colour = flag)) + geom_point()
## Warning: Removed 8863 rows containing missing values (geom_point).

plotly::ggplotly()
all <- rbind(low, hi) |> 
  left_join(wind_dir) |> 
  filter(flag == F)
## Joining, by = "TIMESTAMP"
all <- weatherdash::cat_wind(all, 'wind_dir', 'wind_dir_bin')

ggplot(all, aes(TIMESTAMP, wind_speed, colour = flag)) + 
  geom_point() +
  facet_wrap(~inst_height)

plotly::ggplotly()
#|> filter(as.Date(TIMESTAMP) == as.Date('2022-01-16'))

Calculate \(X_f\):

x_f <- all |> 
  mutate(x_f = -(1/log(f))*(wind_speed/u_star)*((inst_height-d0)/k)) |> 
  filter(x_f < 5000)

Plot \(X_f\) for 3m EC with wind speed on all days: NOTE EC WAS RAISED Feb 17, 2022, this could be why relationship changes after FEB plot.

ggplot(x_f |> filter(inst_height == 3), aes(wind_speed, x_f, colour = wind_dir_bin)) +
  geom_point() +
  facet_wrap(~inst_height+strftime(as.Date(TIMESTAMP), "%m"))

Plot \(X_f\) for 13.5 m EC with wind speed on all days:

Maybe larger footprint during winter months when snow is in the trees?

ggplot(x_f |> filter(inst_height == 13.5), aes(wind_speed, x_f, colour = wind_dir_bin)) +
  geom_point() +
  facet_wrap(~inst_height+strftime(as.Date(TIMESTAMP), "%m"), scales = 'free')

plotly::ggplotly()

Output \(X_f\) table

x_f |> 
  mutate(yr_mon =  strftime(as.Date(TIMESTAMP), "%y-%m")) |> 
  group_by(inst_height, yr_mon) |> 
  summarise(x_f_mean = mean(x_f, na.rm = T),
            u_star_mean = mean(u_star, na.rm = T),
            wind_spd_mean = mean(wind_speed, na.rm = T),
            wind_spd_pk_mean = max(wind_speed, na.rm = T),
            wind_dir_mean = mean(wind_dir, na.rm = T))
## `summarise()` has grouped output by 'inst_height'. You can override using the
## `.groups` argument.
## # A tibble: 14 × 7
## # Groups:   inst_height [2]
##    inst_height yr_mon x_f_mean u_star_mean wind_spd_mean wind_spd_pk_mean
##          <dbl> <chr>     <dbl>       <dbl>         <dbl>            <dbl>
##  1         3   21-11      35.4       0.226         0.613             2.88
##  2         3   21-12      42.8       0.186         0.524             2.85
##  3         3   22-01      55.6       0.186         0.796             2.82
##  4         3   22-02      52.2       0.167         0.629             2.62
##  5         3   22-03      46.4       0.162         0.557             3.74
##  6         3   22-04      53.4       0.141         0.510             1.72
##  7         3   22-05      51.3       0.189         0.588             2.23
##  8        13.5 21-11    1165.        0.500         4.15             10.1 
##  9        13.5 21-12    1118.        0.450         3.53              9.17
## 10        13.5 22-01    1097.        0.524         4.02             11.1 
## 11        13.5 22-02    1061.        0.469         3.39              9.12
## 12        13.5 22-03     957.        0.443         2.96             10.4 
## 13        13.5 22-04     930.        0.420         2.75              9.56
## 14        13.5 22-05    1055.        0.442         3.21              8.71
## # … with 1 more variable: wind_dir_mean <dbl>

Wind Rose for all months for each instrument height:

# for (i in c(3,13.5)) {
#   for (m in 1:12) {
#     fltr <- all |> 
#       mutate(month = as.numeric(strftime(as.Date(TIMESTAMP), "%m"))) |> 
#       filter(inst_height == i,
#              month == m)
#     
#   if(nrow(fltr) != 0){
#     p <- weatherdash::wind_rose(fltr, 'TIMESTAMP', 'wind_speed', 'wind_dir', plot_title = paste('Height', i, "m", 'for', m))
#     file1 <- paste0('figs/wind_rose/wind_rose_height_', i, '_', m,'.html')
#     file2 <- paste0('figs/wind_rose/wind_rose_height_', i, '_', m,'.png')
# 
#     htmlwidgets::saveWidget(plotly::as_widget(p), file1)
#     plotly::save_image(p, file2)
#   }
# 
#   }
#   
# }

Look at one of the windiest days without snow on instruments:

ggplot(x_f |> filter(as.Date(TIMESTAMP) == as.Date('2022-01-16')) , aes(wind_speed, x_f)) +
  geom_point() +
  facet_wrap(~inst_height, scales = 'free')

Wind Rose on this Day:

low_wind <- readRDS('data/eddy_cov_15min_combined.rds') |> 
  filter(as.Date(TIMESTAMP) == as.Date('2022-01-16')) |> 
  select(TIMESTAMP, 
         wnd_dir_compass, 
         u_star, 
         low_ec_rslt_wnd_spd = result_wind_spd) |>
  mutate(
    wind_dir = (wnd_dir_compass + temp_offset) %% 365,
    mean = zoo::rollapply(low_ec_rslt_wnd_spd, width = 48, FUN = mean, na.rm = T, align = 'center', partial = T),
    sd = zoo::rollapply(low_ec_rslt_wnd_spd, width = 48, FUN = sd, na.rm = T, align = 'center', partial = T),
    stdep = (low_ec_rslt_wnd_spd - mean) / sd,
    flag = ifelse(low_ec_rslt_wnd_spd > 4, T, F)
  ) 

ggplot(low_wind, aes(TIMESTAMP, low_ec_rslt_wnd_spd, colour = flag)) + 
  geom_point() 

weatherdash::wind_rose(low_wind, 'TIMESTAMP', 'low_ec_rslt_wnd_spd', 'wind_dir', plot_title = 'test')